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Abstract 

I show stable, localized, single and multi-spot patterns of three classes - stationary, moving, and rotating - 
that exist within a limited range of parameter values in the two-dimensional Gray-Scott reaction-diffusion model 
with cr = D u /D v = 2. These patterns exist in domains of any size, and appear to derive their stability from a 
constructive reinforcement effect of the standing waves that surround any feature. There are several common 
elements - including a spot that behaves as a quasiparticle, a U-shaped stripe, and a ring or annulus, or a portion 
thereof - which combine to form a great variety of stable structures. These patterns interact with each other in 
a variety of ways. There are similarities to other reaction-diffusion systems and to physical experiments; I offer 
several suggestions for further research. 
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1 Introduction 


The Gray-Scott model [2] is a widely-studied model of a pair of reactions involving cubic autocatal¬ 
ysis. It has been applied in reaction-diffusion models in one mmmmm , two pH HQ HE -2Zi |2Sj 
and three [ 2 HI 25 dimensions. Widely-known results include the existence of stable single spots, 
self-replication of spots, spontaneous formation of stripes and hexagonal arrays of spots, and Tur¬ 
ing patterns [9]. Findings that show a lack of moving stable patterns, e.g. [TO] apply only in 1-D or 
a limited region of the 2-D system parameter space ([T6J p. 81 and [19] p. 3). [19] shows existence 
of stable multi-spot patterns within a finite domain. 

The model equations are: 


— = D u V 2 u - uv 2 + F( 1 - u), 

dv o 0 

— = D v V 2 v + uv 2 — (F + k)v. 


( 1 ) 


This paper concerns the two-dimensional case with periodic boundary conditions. I use similar 
terminology and symbols to those in [5]: u and v are the concentrations of two reactants U and V, 
normalized to dimensionless units. Parameters F and k represent the feed rate and removal rate 
of the reactants in the original homogeneous continuously-stirred tank reactor model; in diffusion 
systems they are typically the rates at which U and V permeate through a membrane separating 
a homogeneous supply from the gel in which the patterned reactions occur |6j. D u and D v are 
the diffusion rates of the two reactants through the gel or other medium; considered together as a 
ratio a = D u /D v , they constitute a third parameter that determines certain characteristics of the 
parameter space [9]. I 11 this paper, D u = 2 x 10 -5 and D v = 10~ 5 . Since I am measuring the size, 
velocity, and other related statistics of certain features in observed patterns, there is a need for 
time and length units. For this purpose I use the dimensionless units of length and time implicit 
in the equations, and refer to them as lu (“length unit”) and tu (“time unit”) respectively. 
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2 Methods 


I used numerical simulation similar to that of Pearson [5] (discrete Euler forward integration, also 
called forward-time centered-space m- All figures shown here were produced with a 256 x 256 grid, 
with Ax = 1/143 lu representing a system size of 1.79 lu x 1.79 lu, periodic boundary conditions, 
and At = 1/2 tu. 

Results were verified using higher grid resolution, smaller time steps, and 2 nd and 4 th order 
Runge-Kutta integration. Because many of these patterns exhibit novel properties (such as stability 
combined with motion) it was necessary to perform a thorough series of measurements to isolate 
the legitimate “real” phenomena from any effects that arise from simulation error. 

In these tests, the grid resolution Ax, time step At, and stability index At/(Ax) 2 (the variable 
component of the Courant condition, see M) were each decreased in steps of /2, 2/2 and a/2 
respectively. The finest resolution used was Ax = 1/572 lu, At = 1/128 tu with double-precision 
arithmetic. 


measured velocities range of stable U pattern 


Ax 

At 

(Ax) 2 /At 

U 

3-dot 

(to within ±10 7 ) 

1/143 

1/2 

9.78 x 10 -b 

1 lu/62277 tu 

1 lu/9.1211 x 10 B tu 

0.0608833 < k < 0.0609829 

1/143/2 

1/4/2 

1.38 x 10" 4 

1 1u/62204 tu 

1 lu/9.0037 x 10 6 tu 

0.0608796 <k< 0.0609831 

1/286 

1/16 

1.96 x 10" 4 

1 lu/62300 tu 

1 lu/8.8740 x 10 6 tu 

0.0608778 < k < 0.0609832 

1/286/2 

1/32/2 

2.77 x 10" 4 

1 lu / 62161 tu 

d.n.t. 

d.n.t. 

1/572 

1/128 

3.91 x 10" 4 

1 1 u/62012 tu 

d.n.t. 

d.n.t. 


Table 1: Measurements of the velocities of the two moving patterns in figure 4a performed with progressively 


increasing grid and time resolutions. All tests were performed in IEEE double precision; F = 0.06, k = 0.0609, and 
gridsize 0.9 lu x 0.9 lu or larger. 

Procedure for velocities: 1. Wait for moving pattern to assume its steady-state appearance. 2. Locate a grid 
point on the centerline of the moving feature and on its leading edge, where du/dt is high and it is close to 0.42; 
note its location and the elapsed system time (in lu and tu respectively). 3. Wait for the pattern to travel at least 
200Ax. 4. Measure again, compute v = d/t. Measurement error is 1 part in 200. 

Procedure for last column: 1. F = 0.06; k set to trial value; grid contains a single U pattern. 2. At periodic 
intervals, compute sum of all u values in the grid. 3. Examine the u sums to identify its asymptotic behavior: 
a stable moving U approaches a finite asymptote and continues moving; unstable patterns stop moving or grow 
without bound. 4. Repeat for a new k value; continue until the range of valid k values is known to the desired 
precision. 

d.n.t. = Did not test. 


measured velocities range of stable U pattern 


Ax 

At 

(Ax) 2 /At 

U 

3-dot 

(to within ±10 7 ) 

1/143 

1/2 

9.78 x 10 -b 

1 lu/62420 tu 

1 lu/9.1911 x 10 B tu 

0.0608833 < k < 0.0609829 

1/143/2 

1/4/2 

1.38 x 10" 4 

1 1 u/62477 tu 

1 lu/9.0987 x 10 6 tu 

0.0608796 <k< 0.0609831 

1/286 

1/16 

1.96 x 10" 4 

1 1u/62239 tu 

1 lu/1.0722 x 10 7 tu 

0.0608777 < k < 0.0609831 

1/286/2 

1/32/2 

2.77 x 10" 4 

1 lu/62744 tu 

d.n.t. 

d.n.t. 

1/572 

1/128 

3.91 x 10" 4 

1 lu/62571 tu 

d.n.t. 

d.n.t. 


Table 2: The same tests as in Table 1, performed in IEEE single precision floating-point. 


These tests showed no qualitative difference in the results, but did enable more accurate deter¬ 
mination of certain measurements. For example, the “measured velocities” in Tables 1 and 2 have 
a measurement error of 1 part in 200 due to the procedure used, and the increases in Ax and At 
resolution had no significant effect on the measured velocities: but the double-precision tests show 
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(a) start (b) t = 32 tu (c) t = 120 tu (d) t = 480 tu 


Figure 1: Starting pattern used to generate figure 3c 
and v between 0 and 1, and the background is (u,v) 
exaggerated contrast makes many u values completely black or white. 


In the first image, each rectangle has a randomly chosen u 
(0.4201,0.2878). Coloring is the same as in figure [3j The 


a somewhat greater consistency (i.e. lower variance) in the same measurement across different grid 
resolutions. 

All phenomena reported in the remainder of this paper appear at the grid resolutions of [5], 
Aa; = 1/102 lu, At = 1 tu, with single-precision arithmetic. A few grid effects are apparent at this 


lower resolution, for example the two-spot pattern in figure 4a reorients to a diagonal alignment. 
Such effects diminish sharply at higher grid resolutions. 

The initial state was created in several ways. Most results were produced by starting with 
a background level of u and v set to a homogeneous state computed from F and k. There is 
always a stable state (uhi = = 0), sometimes called the “red state” [5; E, 9]. For (F,k) 

sufficiently small there are two other homogeneous states, one of which can be stable (called the 
“blue state”, due to the color of the pH indicator bromothymol blue the laboratory experiments of 
0)- Transforming the variables and units from Muratov and Osipov [16], the third homogeneous 
state (iih 3 , Vh-s) exists when k < (\f~F — 2F)/2, with u and v given by 


Uh3 — 


A - \M 2 - 4 
2 A 


Vh3 = 


VF(A + VA 2 - 4) 


( 2 ) 


where A is \/F/(F + k). Starting with a grid filled with these values of u and v or with (u — 0, v — 
1), a number of rectangles were added whose width, height, locations, and number were determined 
randomly, and then filled with random levels of u and v. For most tests, these rectangles ranged 
from 5 x 1(T 3 lu 2 to 4 x 10 2 lu 2 in area, and the frequency of occurrence of any given size rectangle 
was inversely proportional to its area. Density of these rectangles ranged anywhere from 1 to 40 
per lu 2 . A typical example is shown in figure [fj 

A repeatable pseudo-random algorithm was used to generate all starting patterns, and the 
initial seed values saved so the same simulation sequence could be reproduced at will m- The 
simulation was run for anywhere from 10 5 to 5 x 10 8 tu, as needed for the phenomena under test. 
A complete survey of the parameter space in figure [2] was performed at intervals of 0.02 in k and 
0.04 in F] many starting patterns were tried at each parameter value combination. The interested 
reader can find a gallery of images online at m- 

Although all of the patterns shown here arose naturally from such initial random states, a 
number of techniques were used to make exploration and discovery more practical: Selectively 
removing unwanted patterns by setting portions of the grid to the homogeneous state; combining 
parts of patterns to create others; changing parameter values; moving small parts of the grid to 
test pattern integrity; then continuing the simulation after making any of these changes. 

Because of their stability, each type 7r pattern naturally evolves from any starting pattern that 
superficially resembles it, provided the dimensions and levels of u and v are approximately correct. 

I found it exceptionally useful to visualize u and du/dt simultaneously via appropriate color 


mappings, with the latter greatly amplified as in figure 4b An interactive simulation tool was 
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essential for the discovery of stable moving patterns “in the wild”. Most work was performed on 
an 8-core Xeon workstation; verification tests were also performed on a PowerPC workstation. 
The simulation was implemented in C, divided into anywhere from 2 to 16 execution threads via 
pthreads. The data set was partitioned in one dimension only (stripes, each overlapping by two 
rows with the two neighbor stripes); all inter-thread data exchange was via shared memory. 

Additional insight was gained from 1-D and 3-D numerical simulations (the latter from the 
website at [21]). 

3 Principal Findings 

I noted all of the pattern types reported in Pearson [5], including several parameter values at 
which two or more of Pearson’s types coexist; many examples are exhibited at CHI- In addition I 
observed three new pattern types, illustrated in figure [3] and named u, £ and 7r to extend Pearson’s 
classification letters. 

Type v is found throughout a large part of the area labeled R in [5], bordering on the regions 
e, A, and //. Here we have stable spots (called solitons or autosolitons by others pT6, [26]) that 
do not multiply. These spots are static only in isolation. Pairs and groups drift apart from one 
another, at a rate that diminishes sharply with distance. The velocity of two spots drifting apart 
in a sufficiently large and otherwise empty domain is modeled fairly well by an exponentially 
decaying rate: velocity c ~ Ke~ rd lu tu -1 where d is inter-spot distance, K and r constant. For 
(F = 0.04, k = 0.07) as shown here, r ^ 46 IlC 1 : the spots drift about half as fast with every added 
0.015 lu of inter-spot distance, r increases with k and with F\ for example, at (F = 0.04, k = 0.072), 
r « 47; and at (F = 0.08, k = 0.066), r ps 67. Given enough time and a small domain size, motion 
ceases with spots roughly equidistant; such multi-spot patterns are the subject of [T9] . At lower 
F values, particularly for F < 0.04, spot-like starting patterns produce spots that oscillate with 
a characteristic frequency varying slightly with k and F, and a damping rate that diminishes 
sharply as (F, k ) approaches the right edge of the type v region in figure [2] Near this frontier the 
stable spot becomes smaller in diameter and the central peak has (u, v) farther from (1, 0). Also in 
figure [2] we see that at higher F values type v occurs on both sides of the saddle-node bifurcation 
line. As Mazin et. ah found in one dimension [9j, at higher F values the red state overpowers the 
blue state in a substantial band to the left of the saddle-node line; however it is not strong enough 
to extinguish small spots. For this reason, parameters in this band support solitons that do not 
grow to fill the space with the blue state. 

Type £ is found in a modest-sized area of the parameter space which includes (F = 0.014, k = 
0.047) and (F = 0.008, k = 0.033). As shown in figure |2j type £ appears on both sides of the 
saddle-node bifurcation. It also appears on both sides of the subcritical/supercritical boundary 
at k — 9/256 ~ 0.035 (when k < 9/256, the model without diffusion oscillates indefinitely; see 
H and [9]). Type £ patterns very closely resemble the B-Z (Belousov-Zhabotinsky) reaction in a 
Petri dish. An initial pattern containing a traveling wave front with free ends produces two spiral 
seeds and usually results in sustained activity. At some parameter values the waves develop flaws 
similar to the segmented waves shown in [23] (which is discussing a CDIMA system); these flaws 
contribute to the production of more spiral seeds. Occasionally, double spirals also occur. A self- 
sustaining population of spiral seeds is needed to maintain the pattern, and the density of seeds 
varies with k and F. The longevity of the pattern is highly dependent on the size of the domain. 
At (F = 0.014, k = 0.047), twenty starting patterns like that in figure [I] (a 1.79 lu x 1.79 lu domain) 
had an average life-span of 5400 tu before all waves died out. Twenty trials in a 6.7 lu x 6.7 lu 
domain with starting patterns of the same density all lasted longer than 10 6 tu. 

Type 7T is the most novel and the subject of the rest of this paper. I found patterns in this 
category for F values ranging from 0.04 to 0.09, in a very thin band of k values; figure [2] shows 
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Figure 2: Parameter values that produce pattern types v, £ and 7 r. The solid curve is the saddle-node bifurcation, 
the dashed curve is the Hopf bifurcation, and the dotted curve is the Turing instability threshold. The area around 
(F = 0.06, k = 0.0609) is indicated by [ 77 ]. 


(a) F = 0.040, k = 0.070 



(c) F = 0.06, k = 0.0609 



(b) F = 0.014, k = 0.047 



u < 0.40 

A 

0.40 < u < 0.44 

■ 

0.44 < u 


(d) shading key 


Figure 3: Pattern types v (inert autosolitons), £ (B-Z targets and spirals), and 7r (stable, stationary and moving 
localized patterns). This color key applies to all subsequent figures. 
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four representative locations in this band. The band runs roughly parallel to the saddle-node and 
Turing lines, and to the left (lower k values) of both. It is also at lower k values than the area 
investigated by McGough and Riley EE]. 

In [5] Pearson explored some parameter values near this area, initializing the system with red 
state (u — 1,d = 0) and a small central rectangle of (u = 1/2 ,v = 1/4), and found that the 
system evolved to a homogeneous blue state; these were designated B. Mazin et. al. |9j using 
a similar starting pattern in one-dimensional simulations, explored a greater number of F and k 
combinations, and found many for which stable localized structures result; these were designated 
L 

At all F values the n band is very narrow in k as compared to its distance from the bifurcation 
lines. For example, when F = 0.06, the saddle-node bifurcation is at k ~ 0.06247, the Hopf 
bifurcation is at k ~ 0.06245, the Turing bifurcation is at k ta 0.06191, and the range of valid k 
values for the object in the lower-left of figure |4a| is 0.06087 ± 0.00001 <k< 0.06098 ± 0.00001. 
Thus the width of the 7r band near F = 0.06 is about 1.1 x 10~ 4 in units of k, while the distance 
to the nearby Turing region is ten times as great. 

Pearson’s pattern type t (see figure 2 in m shares some qualities with type 7r. It is in region 
L of Mazin et. al. [9], but its lack of solitary spots indicates it is probably slightly outside the 7r 
band. 


4 Stability and motion of type 7r Patterns 


Figure [4a] shows five patterns at F — 0.06 and k = 0.0609. The three in the upper row are stable 
non-moving patterns. Contrast has been exaggerated in the area of 0.40 < u < 0.44 to show the 
concentric rings or “halos” that surround all patterns in systems with these parameters. The halos 
are concentric stationary waves of alternating sign superimposed on the homogeneous state values 
[uh'ii Vh'i) in (j2]). The full range of u and v levels for type n patterns at F — 0.06, A: = 0.0609 is 
0.29 < u < 0.86 and 0.01 < v < 0.43. 


The other two patterns in figure 4a move to the left, indefinitely at constant speed; the three- 


spot pattern moves at about 1 lu per 8.5 x 10° tu, and the U-shaped pattern at about 1 lu per 
6.2 x 10 4 tu. 

All of the patterns in [4a] arise frequently from random starting patterns, and are resilient to 
noise and other perturbations. If any of them is perturbed by shifting half of the pattern in any 
direction a distance on the order of 0.02 lu, further simulation results in a return to the canonical 
forms shown here. 

Figure 4b shows the time derivative of u for the same five patterns. The motion of the two 
patterns in the bottom half of the figure is clear. 


Figure 4c shows two rotating patterns. Both rotate clockwise; the four-spot pattern performs 


one full revolution in about 1.6 x 10 7 tu; the other takes about 1.2 x 10' tu. Disruptions to these 
will cause a momentary change in rotational velocity followed by a return to the normal rotation 
rate after the spots return to the stable alignment. 

Many stripelike patterns like those in figure [4d] arise from random starting patterns; one-ended 
forms are more common than the two-ended versions shown here. They are stable in the central 
linear section, in the direction perpendicular to their length, but are unstable in the other dimen¬ 
sion. Any deviation from a straight line will increase, first slowly and then with increasing speed, 
forming meanders like those of a river in a floodplain. These linear forms grow at both ends. The 
top one is the fastest, each end grows at about 1 lu per 6.7 x 10 4 tu, slightly slower than the 


speed of the U-shaped pattern in figure 4a The second example grows at 1 lu per 1.45 x 10° tu 


from each end. The third example has two spots at each end that make it grow (at about 1 lu 
per 2 x 10 6 tu); without these spots the pattern shrinks at each end. When encountering other 
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(a) u (b) du/dt (c) rotating patterns (d) stripelike patterns 


Figure 4: (a),(b): Comparison of u to its time derivative for pattern type 7 r. Coloring for all but (b) is the same 
as above. In image (b), white represents du/dt > 2 x 10 -6 , shades of gray where \du/dt\ < 2 x 10~ 6 , and black 
otherwise, (c): Two stable rotating patterns, (d): Three linear growing patterns. (F = 0.06, k = 0.0609) for all 
figures; details in section |4j 


objects, these stripelike patterns will change direction, stop growing, or (frequently in the case of 
the topmost example) break down into some other form, such as separate parallel stripes. 

Systems that contain active growing stripes will usually grow to fill all available space with 
spots and stripes, but the time taken and the final proportion of spots to stripes is highly variable, 
and depends on tiny details of the initial configuration. This is substantially different from the 
behavior of stripes at other parameter values in the Gray-Scott system, whose final density can 
generally be predicted from the parameter values alone. 

The central portions of the stripelike patterns in figure [4d] suggest a connection with certain 
one-dimensional localized structures found by Mazin et. al in |J] at parameter values in their region 
L. In my tests 1 found that the dimensions and levels of u and v of features in 1-D simulations were 
effectively the same as cross-sections of stripe-like patterns in 2-dimensional simulations when the 
same parameter values are used. 


5 Variations with small changes to k 


The patterns in figure [4a] are shown as they appear when F = 0.06 and k = 0.0609. If k is 
diminished below 0.06087, the U-shaped moving pattern is no longer viable and decays into a 
single spot. When k < 0.06062 the single spot shown in the upper-left of figure 4a is no longer 
stable and quickly evolves to the homogeneous state, however the other three patterns made up of 
spots continue to exist. When k is lowered below 0.06060, the patterns in the top center and lower 
right vanish, but the triangular three-spot pattern remains. Below 0.06057 the triangular pattern 


vanishes, but a 7-spot arrangement similar to that in figure 5c is still viable. The 7-spot pattern 
dies out when k < 0.06055. 


At F — 0.06, k = 0.06090 the 5-spot pattern in figure 5a is stable and moves to the left slowly. 
Using this as the starting pattern, k is increased to 0.06110, and a spot forms in the center to 
make the 6-spot pattern in figure 5b. This pattern is also left-moving and stable. Increasing k 


further to 0.06135, a seventh spot appears (figure 5c), and motion stops due to the attainment of 
symmetry. Throughout this process the spots shift away from one another slightly to attain a new 
steady state, with a slightly greater inter-spot spacing with each increase in k. 

Increasing k to 0.06150, new spots appear on all sides, quickly growing to produce the uniform 
hexagonal grid in figure [5d| 

Groups of spots are no longer stable at these relatively high k values — if A; is increased more 
gradually from 0.06135 to 0.06150 (for example in 15 equal steps at intervals of 5000 tu), some or 
all of the spots will have time to swell into elongated “stripes” of high u. Even at k — 0.06135 


the seven-spot pattern shown in 5c is barely stable, and groupings of more than 7 spots produce 
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(d) k = 0.0615 


(e) k = 0.0620 


Figure 5: (a) through (d) show the evolution of a starting pattern as k is raised progressively to three nearby values. 
For comparison, (e) shows a Turing pattern arising from low-level random noise. F = 0.06 for all figures; details in 
section [5] 


stripelike patterns. 

Considered together, the range of different behaviors described thus far represents a very wide 
spectrum of phenomena in the relatively narrow range 0.06055 < k < 0.06150. A similar spectrum 
is found at higher and lower F values, although varying in details. For example, the behavior 
shown in figure 5d has been observed at several values of F from 0.046 to 0.0652; at higher F 
values a greater change in k is needed to precipitate the growth of the hexagonal pattern, and it 
appears that at significantly higher F the phenomenon cannot be produced at all. The U-shaped 
moving pattern has been found to be stable at F values from 0.0492 to 0.0876, with k varying as 
illustrated in figure |2j 

At F — 0.06, k = 0.062, the Turing effect is present. Here, no initial spots are needed to generate 
a pattern of spots and/or stripes. As defined by [I] (see also the introductory section of |9j) Turing 
patterns arise spontaneously from random noise of arbitrarily low initial amplitude, diffusion plays 
an active role in the destabilization of the initial state, and a characteristic wavelength (spot size 
and/or stripe width) exists that is independent of the system size. Figure |5e] shows a typical Turing 
pattern in the 2-D Gray-Scott system; it was produced from a starting pattern of “white” noise 
with amplitude 1CT 4 superimposed on the blue state (■ Uh 3 ,Vh 3 ) from (|2]) . We see similarities, such 
as hexagonal arrangement of spots and a similar spot size and spacing, to figure 5d However, 
many grain boundaries, short stripes and other differences in detail are present. This is typical 
of patterns throughout the Turing domain, which at F = 0.06 includes k values in the range 
0.06191 < k < 0.06245 (see pj for a derivation of the formulas for these values; k = 0.06245 is the 
Hopf bifurcation). 


6 Distance interactions and more exotic patterns 

Figures [6a] and [6b] show two frames of a simulation in which k was gradually increased from 0.06058 
to 0.06110 over 10 6 tu. The overall length of the horizontal row of spots increases with k, and 
the smaller 3-spot pattern moves slightly away from it. The entire set of 15 spots also rotates 
very slowly as a unit. Similar interactions between patterns that are not in direct contact are very 

8 



















0 ©^ 


(c) Interaction (before) 

Figure 6: (a),(b): Lattice spacing variation; (c),(d): An interaction between three patterns. (F = 0.06, k = 0.0609) 
except as otherwise noted; details in section [6] 



Figure 7: (a)-(c) are precursors to (d)-(f), which are stable and rotate clockwise. Patterns in (g) move to the right; 
the asymmetrical pattern moves on a curved path. Details in section |6l (F = 0.06, k = 0.0609) for all figures. 



(d) Interaction (after) 



common. Decreasing k back to 0.06058 over a similar time period causes the pattern to return to 
the state inIHal 

In figure 6c, the U-shaped feature is moving to the right and slightly downward. It interacts 
weakly with the stripelike feature, then more strongly with the triangular array of spots to the 
right, and emerges traveling on a diagonal path towards the upper-right. An animation of this is 
available at [ 28] , 

In figure [7] are six rotating patterns with a “target” pattern in the center. The first three are 
asymmetrical; each of these transforms into the corresponding symmetrical version; the center 
shifts a bit in the process. All rotate clockwise. The solitary spots in the second example are 
pushed by the triads. In figure [7b] one of the spots is stationary; after becoming symmetrical as 
shown in 7e both spots are moving. Each of these solitary spots remains at the distance shown, 
about twice the radius of the brightest surrounding halo, from the annulus of the “target” and 
from the two closer spots of the triad immediately following it. 

Figure 7g shows several more moving patterns. All move to the right. The pattern in the 
upper-left is an example of many variations of the “target” with attached spots that arise from 
random starting patterns. Most remain stable and move; symmetrical forms move on a straight 
path and asymmetrical forms move on a curved path. The “target” by itself is also stable and does 
not move. The other three all involve spots that are being “pushed” ahead of a larger pattern; 
in all cases the larger pattern moves faster if the spots in front are removed. The asymmetrical 
pattern in the lower-right moves on a curved path, bending to its right as it moves forward. In the 
process the spot is pushed out of the way, and is left behind after about the first quarter-turn of 
the moving pattern’s path. An animation of this is available at 
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7 Similarity to Other Systems 


Several other authors describe spots with appearance and interaction similar to the spots shown 
in this paper. 

Schenk, et. ah [12], studying a two-component cubic autocatalytic reaction-diffusion system, 
found parameter values for which there is a stable spot solution with concentric rings of alternating 
sign and diminishing magnitude. For certain parameter values they found that the spots interact 
in a similar way to that described in this paper; they showed some stable multi-spot arrangements 


similar to my figure [4a] and proposed several more. There are some differences, for example the 
pattern in their figure 5(e) (a four-spot pattern in a Y configuration with unequal angles) is not 
stable in the Gray-Scott system with parameters (F = 0.06, k = 0.0609), nor at several other 
locations I tested throughout the n region. I also found the pattern in their figure 5(b) (three 
spots in a straight line) to be stable, whereas they did not. Such differences might result from 
differences in the spacing and relative amplitudes of the halos. 

In work related to [12], a number of researchers including Purwins et. al. [13] have identified 
and studied “quasi-particles” that can arrange themselves in a very similar manner to my figures 
This can be seen clearly in figures 3 (d), (g), and (h) of the web page article by 
Note in particular the halos around the spots in their figure 3(d). These spots 


4a and 5d 


Stollenwerk 

arose in laboratory AC gas plasma discharge experiments that have other features reminiscent of 
patterns seen in the two-dimensional Gray-Scott system. 

Liehr, et. al. in [T5] . modeling a DC gas plasma discharge experiment, began with a 3- 
component reaction-diffusion system, then transformed the equations into a two-component form 
with a global integration term. In numerical simulation (see their figures 3 and 5) they achieved 
results similar to my figures 5c and 5d Note the halos and arrangement of the spots in their 
patterns. Studying a similar 3-component system, Schenk in [44] shows several examples that 
include a Zielscheiben-Struktur (target structure) closely resembling the target patterns in my 
figure [7} Spots with halos and similar spot-to-spot interactions are also shown in that work. 
In both of these 3-component systems, a single spot can have an intrinsic velocity and display 
particle-like collision behavior in addition to the static attraction and repulsion effects. 


8 Discussion 

The single spot in the upper-left of figure [4a] has a set of concentric halos with progressively lower 
amplitude and alternating sign. In multi-spot patterns, each spot tends to be found at a location 
that coincides with the first positive-sign halos of neighboring spots. In the growth of hexagonal 
arrays, the new spots always appear at such locations. 

In figure [44] the U-shaped moving pattern is of similar size and shape to the three-spot pattern, 
and both tend to return to the canonical dimensions shown here after distortion or perturbation. 

In more distant alignments such as that in figure [6a] the spots are found at distances coinciding 
with each other’s second positive-signed surrounding halos. This is also evident in the five-spot 
pattern in figure [5a] which remains as shown when perturbed, rather than changing to a pentagonal 
ring or some other arrangement. 

In several instances (for example figure 4d, all of figure [6] and figure [7d| we see rows of spots 
and other roughly linear features, halos that parallel these features, and motion and interaction 
that follows the locations of these roughly linear halos. 

All of the foregoing suggest that areas of high u and low v produce a pattern of surrounding 
standing waves, and that these standing waves combine in a nearly linear way to produce an effect 
of constructive and destructive interference, causing spots and other features to drift towards a 
preferred alignment. A similar explanation exists for spot interactions in the reaction-diffusion 
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system of ra- Analytical research is needed to establish a basis for this theory for the Gray-Scott 
model equations, or another explanation that can account for the observed phenomena. Inasmuch 
as similar effects are observed in the 1-dimensional system (see Mazin et. al [9] and my own work 
m, analytical work can probably begin with 1-dimensional systems. 

There is much opportunity for further research by numerical simulation, including rigorous 
statistical analysis to establish the degree of stability of the patterns in response to varying types 
and levels of perturbation. The majority of the present results are at a single point in the parameter 
space; there is a likelihood that changing F and a = D u /D v will yield new discoveries. 

The superficial similarity (spot shape and “halos” and target patterns) and behavior (spots con¬ 
densing into a bound multi-spot state similar to “molecules”) in gas plasma experiments, combined 
with the successful modeling (e.g. by [15]) of these experiments by a reaction-diffusion system in 
numerical simulation, strongly suggest research to identify more connections between the present 
work and the gas discharge systems. 

Concerning the moving stable structures, there are several well-studied reaction-diffusion system 
models, including the Brusselator model with diffusion, the three-component system of [15], and 
the Oregonator model of the B-Z reaction |18j : all sharing many or most of the features of the Gray- 
Scott system (stripes, spots, multiple homogeneous stable states, mixed-mode patterns, ability to 
produce Turing patterns, and varying types of movement and interaction at various parameter 
values). Because of these many shared traits, it seems likely that some of these other systems can 
also produce stable moving localized patterns and other complex phenomena like those described 
in the present paper. If such phenomena exist, they likely involve parameter values that have not 
been thoroughly explored, or may be confined to a very narrow region of the parameter space. 

Finally, this author notes that the great diversity of patterns and types of interaction displayed 
by the Gray-Scott system at these specific parameter values, combined with the inherent stability 
of these patterns, clearly places this system in Wolfram’s class 4 (complex localized structures, 
sometimes long-lived) [3] that characterizes certain discrete cellular automata that have been 
shown to be capable of universal computation. Constructing large systems of interacting patterns 
presents several challenges for long-term stability, because most interactions cause each interacting 
part to shift. However, a starting pattern containing an infinite array of elements, each of which 
is only used a small number of times, would appear promising. 
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